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■ A restricted planar circular three-body system, consisting of the Sun and two planets, is 
studied as a simple model for a planetary system. The mass of the inner planet is considered 



to be larger and the system is assumed to be moving in a uniform interplanetary medium 



OO - 

Q\ . with constant density Numerical integrations of this system indicate a resonance capture 



when the dynamical friction of the interplanetary medium is taken into account. As a 
result of this resonance trapping, the ratio of orbital periods of the two planets becomes 



O 

nearly commensurate and the eccentricity and semimajor axis of the orbit of the outer 



planet and also its angular momentum and total energy become constant. It appears from 
the numerical work that the resulting commensurability and also the resonant values of 
the orbital elements of the outer planet are essentially independent of the initial relative 
positions of the two bodies. The results of numerical integrations of this system are 
presented and the first-order partially averaged equations are studied in order to elucidate 
the behavior of the system while captured in resonance. 
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1 INTRODUCTION 

The dragging effect of the interplanetary medium and its role in formation and evolu- 
tion of planetary systems are quite well known and have long been investigated by many 
authors such as Weidenschilling et al. (1985 k 1993) , Patterson (1987) , Peale (1993) , 
Malhotra (1993) , Beauge et al. (1993 and 1994a&b) , Murray (1994) and Gomes (1995). 
Recently, it has also been shown that the dynamical friction of this medium can account 
for the stability of planetary orbits in the early stages of their dynamical evolution (Melita 
and Woolfson 1996, here after MW). Using a formulation developed by Dodd and McCrea 
(1952, see also Binney and Tremaine 1987) and by integrating a general three-body system, 
MW showed that the system of Sun and two massive planets, subject to dynamical friction 
and accretion, will be captured into a near mean-motion resonance when the inner body is 
more massive. As a result of this resonance trapping, the eccentricities of the both planets 
become constant and their semimajor axes decrease monotonically. 

A simple planetary model is studied here based upon the results of Melita and Woolf- 
son (1996). The simplifications introduced in this study make it possible to analyze the 
averaged dynamics of the system, analytically. The model used in this paper is a restricted 
planar circular three-body system with an inner planet that is more massive. Since MW 
showed that the effect of accretion will not change the end results qualitatively, this effect 
is not included in the dynamical equations of the system here. However, the frictional 
effect of the interplanetary medium is the main source of non-gravitational dissipation and 
is therefore taken into account. 

The model under investigation is presented in section 2. Section 3 is concerned with the 
results of the numerical integrations. In order to analyze the numerical results analytically, 
a newly developed averaging technique for dissipative systems is introduced in section 4. 
This averaging method has been developed by Chicone et al. (1996-7a&b) and was used 
in their analysis of dynamical behavior of binary systems subject to external gravitational 
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radiation as well as radiation reaction damping. This averaging technique provides a very 
helpful tool to study the dynamics of the system at resonance. In section 5, the application 
of this method to the dynamical equations of the system are discussed. The first-order 
averaged system at resonance is presented in section 6 and in section 7, conclusions of this 
study are presented by reviewing the results and discussing the possible extensions and 
applications. 



2 THE MODEL 

A restricted planar circular three-body system consisting of a central star (hereafter 
S) and two planets is considered with the inner planet more massive. It is assumed that 
the gravitational attraction of the planets on the central star is so small that its motion 
is negligible. This allows for the consideration of an inertial coordinate system with its 
origin on S . It is also assumed that after taking all perturbative effects into account, 
the orbital motion of the inner planet (hereafter X) is uniformly circular with a given 
period. Dynamics of the outer planet (hereafter O) is thus the focus of attention and is 
determined by the gravitational attractions of the central star and the inner planet and 
also by dissipation due to dynamical friction caused by the interplanetary medium. 

In the inertial coordinate system under consideration, the equation of motion for 
planet O can be written as 

d2 ? r Mm _ m I m ^ 

m o —rpr = ~Q — 3— r Q - Q —, ^ r D - rj + R f , 1 

at 1 r A Q \r Q — r T \ 6 

where the indices I and O stand for the inner and the outer planets respectively, and ^ / 
represents the dynamical friction force of the interplanetary medium and is given by (Dodd 
& McCrea 1952, Binney & Tremaine 1987 ) 

^--^^(.♦g)* . (2) 
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In this equation, W is the relative velocity of O with respect to the medium, p is the 
uniform density of the medium and 

S = r , 0) 

where M represents the mass of the central star. 

Equation (1) can be simplified by dividing both sides by m and also by letting 
^ = T$f/m and r = r Q . The main equation under consideration can therefore be 
written as 

w = -g-,f-g w -^(f-f,) + ^ . (4) 

It is more convenient to write this equation in dimensionless form. Introducing t and r 
as the quantities that carry units of time and length respectively, t and r can be written 
as t = t t and r = r r , where t and f are their corresponding dimensionless variables. 
Equation (4) has now the form 

* ? = -K*-eK^U + l , (5) 



dt 2 ' f 3 |f - f,! 3 



where 



K = . (6) 



^0 



£ = — (7) 



and 



t 2 



a = 3 e) . (8) 



In the rest of the calculations, I choose a set of units such that K = 1 . Therefore from 
equation (6), i and r will be related as t 2 = rfJQM . This relation implies that r can 
be viewed as the orbital radius of a planet that orbits S on a circular path with a period 
27r£ . In the model presented in this paper, planet X has this type of motion. Therefore, 
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in the rest of this paper, I set r = r T and t = Tjlis where T T represents the (given) 
orbital period of planet X . Dropping the hat signs, equation (5) can be simplified as 



d 2 f f (f — ?j ) 

dt 2 r 3 \f — ?j | 3 



+ 1 , (9) 



where r I is the unit vector along r I . 

Let us now assume that the motion is planar. Introducing a polar coordinate system 
on the plane of the orbits and with its origin at the location of the central star, the 
dynamical equations of planet O will be given by 

P r = r , (10) 



Pe = r 2 6 , (11) 



p r = ~ ^, ~ [f |3 [r - cos (9 -9;)] + (R x cos9 + R y sm9) , (12) 



and 



P e = - e T ^ s\n{9 - fl f ) + r(- R x sing + R y cos 6) , (13) 



\r — r z I" 3 



where R x and R y are the x and y components of Rl , respectively, and 9 X = 2nt/T I + 
constant. In dimensionless form and with a time origin such that this constant value 
becomes zero, we have 9 t = t , where t is the dimensionless time variable. 

In order to calculate R x and R y , we need to turn our attention to the interplanetary 
medium. It is assumed that this medium is freely rotating around S (Kiang 1962 , Dormand 
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& Woolfson 1974). Therefore its velocity at any point is perpendicular to the position 
vector at that point with a dimensionless magnitude equal to = r -1 / 2 . The relative 
velocity W = v — , will therefore have a radial component equal to r and a tangential 
component equal to r9 — ru^ where uj^ = r -3 / 2 is the dimensionless angular frequency 
of the medium at distance r . As a result, the magnitude of W will be equal to 

W 2 = r 2 + r 2 (8 - w M ) 2 . (14) 

The components of W are easily obtained noting that at any distance r , the angle between 
Vfj, and the x-axis is equal to (6 + f) . Therefore 

W x = f cos 6 — r (6 — ojn) sin6> , (15) 

and 

W y = r sin6 + r (6 — u; M ) cos^ . (16) 
From equation (2) and in dimensionless form, R x and R y will then be equal to 

R x = - A- In (1 + Br 2 W 4 ) [fcos6 - r(6-uj^sm6] , (17) 

and 

R y = - A. i n (i + Br 2 W 4 ) [f sin 6 + r(6 - oo^) cos 9] , (18) 
where A and B are positive parameters given by 

such that A « 1 and AB « 1 for any realistic system. 
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3 NUMERICAL RESULTS 



Equations (10) to (13) were numerically integrated for different values of the pertur- 
bation parameter e and density p using a stiff integration routine. Following MW, inte- 
grations were first performed on the Sun- Jupiter-Saturn system starting from the present 
near (5:2) commensurability. The density of the interplanetary medium was taken to be 
constant and equal to 10 -11 Kg m -3 which is equivalent to 16 times the mass of Jupiter 
uniformly spread in a sphere of radius 50 au. In complete agreement with the results of 
MW, a near (2:1) commensurability is obtained for all initial relative positions of the two 
planets. A typical case is presented in Figure 1 . 





2 . b[ a 
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Figure 1. Graph of the ratio of orbital period of the outer planet (i.e., Saturn) to that of the inner one (i.e., Jupiter) 
against time. Integrations were performed with Jupiter on the x-axis with 9 j = and Saturn at 9 Q =45° . The 
initial values of T , Pg and P r were calculated from r = Cl(l — e 2 ) /(l + e COS v) , Pq = a(l — e 2 ) and 
P r = e sin V / Pg , where 0=1.838046 is the present dimensionless value of the mean semimajor axis of Saturn's 
orbit, e = 0.0556 is its present mean orbital eccentricity and V = 15° is its initial true anomaly (see section 4 
and also Figure 4). The timescale is (10 4 T / /27T) years. 
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As a result of this resonance trapping, the orbital eccentricity and semimajor axis of 
the outer planet (i.e., Saturn) and also its angular momentum and total energy become 
essentially constant (Figure 2). 




1 2 3 4~ 




1 2 3 4 5 



Figure 2. From top to bottom, orbital eccentricity and semimajor axis (left column) and angular momentum 
and energy (right column) of Saturn against time while captured in resonance. The initial conditions and the 
timescale are similar to those of Figure 1 . 



Since the present density of the interplanetary medium is much smaller than 10 _11 Kgm -3 , 
integrations were also performed with lower values of density for the actual Sun- Jupiter- 
Saturn system. A typical case is shown in Figure 3 , where the system is captured in a 
(3:1) resonance but leaves this state after a short time. 
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Figure 3. Graphs of the Ratio of orbital periods (left) and orbital eccentricity of Saturn (right) versus time for 

1 A Q 

the actual Sun-Jupiter-Saturn system in an interplanetary medium with a density equal to 10 Kg m 
The system is captured in a (3:1) resonance and leaves the resonance after a short time. The initial conditions of 
integration are given by (a, e, 9, v) = (1.838046 , 0.0556 , 75° , 30°) and the timescale is (10 5 T / /27r) years. 

In another series of runs, the planetary masses in the Sun-Jupiter-Saturn system were 
changed while keeping their mass ratio constant and equal to the present mass ratio of 
Jupiter and Saturn. The density of the interplanetary medium was kept equal to the 
16 times the mass of the inner planet uniformly distributed inside a sphere of radius 50 
au and the mass of the central star was unchanged and equal to that of the Sun. Inte- 
grations were performed for different values of e and it was observed that only for the 
values of e in the range 10~ 2 to 10~ 4 , the system was captured in resonance (Figure 4). 
The interesting result was that starting from present near (5:2) commensurability between 
Saturn and Jupiter, all resonances for this range occurred near (2:1). It is important to 
mention that for the values of e smaller than 10 -4 , the planetary model under considera- 
tion will no longer be physically valid. Numerical results indicate that in these cases, the 
amplitude of oscillations of r grows rapidly with time which results in change of the con- 
figuration of the system. That is, for long amplitudes, planet O approaches X and crosses 
its orbit. This results in an exchange of positions between O and X and also causes the 
quantity \r — r x \ in the dynamical equations of the system, reaches zero at certain times. 
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Figure 4. Graphs of the ratios of orbital periods (left column) and orbital eccentricities (right column) versus 
time, for 6 = 10~ 2 (top), and 6 = 10 -4 (bottom). The initial values are given by (a, e, 6>, v) = (1.838046 , 
0.0556 ,75° , 30°) and the timescale is (10 4 T 7 /27T) years. 



Numerical experiments on the actual Sun-Jupiter-Saturn system, were also performed 
with larger values of the mass of the central star M , but no resonance capture was obtained. 
A typical case is presented in Figure 5. 
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Figure 5. Ratio of orbital periods (left) and orbital eccentricity of the outer planet (right) against time for a 
central star with a mass equal to 1000 times the present mass of the Sun. The initial values and the timescale 
are the same as those in Figure 4. 
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The rest of this paper is devoted to the analysis of the results of the numerical in- 
tegrations presented in this section. For this purpose, I will use a method of averaging 
called "Partial Averaging Near a Resonance" (Chicone et al. 1996) which is described in 
the next section. 

4 AVERAGING 

Consider the following perturbation problem in a /c-dimensional Euclidean space K fc , 

u = /(«) + eb(u,t,e) , u G K fc . (20) 

In this system e is the perturbation parameter, b is a periodic function of time with a 
period r\ > and the overdot represents the derivative with respect to time t . If the 
unperturbed system, i.e. u = f(u) , u G K fc is integrable, one can always find a set of 
canonical coordinate transformations to write equation (20) in terms of a set of action- 
angle variables (£, </>) G 3? fe x W- such that 

t = eT{CA) + 0(e 2 ) , (21) 

and 

= 0(£) + e £(£,</>) + 0(e 2 ) . (22) 

In these equations is a 2tt modulo ^-dimensional vector of angular variables and T and 
B are 2tc periodic functions of <fi . The most important feature of equations (21) and (22) 
is that in the angular equation (22), there appears a term which is independent of the 
perturbation parameter e. This indicates that, in the first order of approximations (i.e., e), 
4>(t) can be considered as a fast-changing angular variable while C(t) moves slowly away 
from its unperturbed value. Therefore, when studying the dynamics of the system (20), 
the time evolution of the averaged value of C(t) over <f>(t) can be taken as a reasonable 
approximation to the evolution of C(t). 
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The above statement is often referred to as the "Averaging Principle." This principle 
states that if C = Co is the initial value for the action variable C in system (21) and (22), 
then for sufficiently small e , the solution to the initial value problem 



will provide a useful approximation for the evolution of the action variable C with an error 
of order e 2 over a timescale of order e _1 . In this equation T is the averaged value of T 
and is given by 



In the system presented in this paper, the unperturbed system (i.e. Sun-Saturn) is 
Hamiltonian. Also, in the perturbation problem given by equations (10) to (13), e and A 
are small quantities. Therefore we can use the averaging principle to analyze the dynamics 
of this system while captured in resonance. To this end, one needs to write the equations 
of motion in terms of appropriate action-angle variables that will provide us with a fast- 
changing angular variable. This is accomplished by writing the equations of motion in 
terms of Delaunay's variables which can be described as follows. 

If at any time t all perturbative forces are removed, the orbit of the outer planet will 
become an ellipse with its focus always at the origin of coordinates. This ellipse is tangent 
to the actual orbit at point r(t) and is therefore called an osculating ellipse (Brouwer & 
Clemence 1961, Kovalevsky 1967, Hagihara 1972 and Gutzwiller 1998). One can use the 
orbital elements of the osculating ellipse to introduce the Delaunay action-angle variables. 
In our planar problem, the relevant Delaunay variables are given by (Kovalevsky 1967, 
Sternberg 1969, Hagihara 1972) 



J = eF(J) 



J(0) = C 



(23) 




(24) 



L = a 1 ' 2 , 



(25) 



I = u — esin-u 



(26) 
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G= [a(l-e 2 )] 



2M 1/2 



(27) 



and 



9 = 



(28) 



where a and e are the semimajor axis and eccentricity of the osculating ellipse, respectively, 
and / is the mean anomaly. In these equations L and G are action variables, while / and 
g are their corresponding angular variables. The variables v and u are, respectively, the 
true anomaly and the eccentric anomaly of the outer planet as illustrated in Figure 6. 

The polar coordinates and the momenta of the outer planet are related to the Delaunay 
variables by Pg = G 2 = r(l + ecosO) and P r = e smv/G . Using these equations along 
with replacing 6 from equation (28), the equations of motion in terms of the Delaunay 
variables can be written as 



dL 
~dt 



a(l-e 2 ) 



-1/2 



F r e sin-0 + Fg (1 + e cost>) 



(29) 



dG 



dl 

— = u + - a 

dt e 



F r [— 2e + cos-0 + ecos 2 v) — Fq (2 + e cos-D) sin{i 



(30) 
(31) 



and 



where 



dg_ = 1 

dt e 



a(l-e 2 ) 



1/2 



-F r cos v + Fq(i H r^sin-O 

V 1 + e cos v / 



F fl = -e 



r — r T 



sin(^-^; 



r — r T 



(i? x sin — R y cos 0) , 



(32) 

(33) 
(34) 



and the Keplerian frequency of the osculating ellipse is given by u> = L 



-3 




Figure 6. The osculating ellipse for the outer planet. The central star is at one of the focal points. The orbit of 
the inner planet is a circle. The true anomaly V and the eccentric anomaly U are defined according to this figure. 



Quantities F r and Fq are, in fact, representing perturbations due to the gravitational 
attraction of the inner planet and also the dynamical friction of the interplanetary medium. 
It is evident that only the angular variable / , i.e. the mean anomaly, has an associated 
Keplerian frequency unaffected by perturbations. This immediately suggests that / is 
the appropriate "fast" angular variable for averaging purposes. In the next section, the 
averaging principle is applied to the system of equations (29) to (32) and the averaged 
dynamics of the system while captured in resonance is discussed. 
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5 AVERAGED SYSTEM AT RESONANCE 

As mentioned above, equations (29) to (32) represent the dynamics of the system in 
terms of action variables (L,G) and angular variables (l,g) ■ The objective of this section 
is to apply the averaging principle to the equations (29) to (32) while the system is at 
resonance. In general, resonance occurs when the Keplerian frequency uj becomes com- 
mensurate with the frequency of the external perturbation. In the system presented in 
this paper, the external perturbation is due to gravitational attraction of the inner planet 
and its frequency is given by Uj = 2 , k/T i . Therefore in this system, resonance capture 
means that relatively prime integers m and n exist such that mu = nuj 1 . Since in this 
model the inner planet is orbiting the central star on a circular path with a fixed orbital 
period, the immediate result of resonance trapping is a constant value for the orbital pe- 
riod of the outer planet. A constant orbital period in turn, results in a constant value for 
the semimajor axis of the osculating ellipse and therefore a constant value for the action 
variable L . At resonance, L = L$ = (m/nu;,) 1 / 3 . However, as it is evident from Figure 2, 
the value of L = a 1 / 2 during resonance is not entirely constant, but oscillates around its 
Keplerian value L . Let D be a measure of these oscillations and also let the oscillations 
of the mean anomaly from its Keplerian value Lq 3 t be given by (p . Requiring L and / to 
have the same power of e in the lowest order of approximation, one can write 

L = L + e 1/2 D , (35) 

and 

1= j^t + cp , (36) 

where the first order of approximation here has perturbation parameter e 1 / 2 . 

Equations (35) and (36) are indeed the necessary transformations that one needs to 
employ in order to write equations (29) to (32) for the system at resonance. In order to 
study the averaged dynamics of the system at this state, the averaging method presented 
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in section 4 must then be applied to the transformed equations. However, it is important 
to mention that the averaging principle is valid only for systems with one angular variable 
(Chicone et al. 1997a). In the system presented in this paper, there are two angular 
variables / and g . One can, however, show that the average rate of variation of g is in fact 
negligible and therefore it can be considered as a slow-changing quantity. This becomes 
evident by writing equations (29) to (32) as 

dL dH ext 

-dF = - £ ^r + eAnL ' (37) 

dG dH PT f A „ ,„„ N 

7F = - £ ^f + eAKc - • (38) 

dl 1 dH PT f . „ ,„„ N 

* =u + £ ^f + eAK ' ' (39) 



and 

dt ' OG 

where A = AB/e , 



do dH ext . ^ . , . 



H ext = - 1 ^. , (41) 
\r — r r | 

and 7?-l , , 7^.^ and lZ g represent the contributions of dynamical friction in their corre- 
sponding equations. I will show in section 6 that from these contributions, only TZl will 
appear in the first-order averaged system. 

Equations (39) and (40) indicate that while / has a perturbation-free frequency L~ 3 , 
the rate of change of g is proportional to e . This implies that, in comparison to / , the rate of 
variation of g can be neglected and therefore, for the purpose of this study, equations (29) to 
(32) can be considered to have only one fast-changing angular variable that is / . Equations 
(37) and (39) have now the general form of the system (21) and (22) and can therefore be 
used as the appropriate grounds for applying our averaging method. That means, after 
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replacing L and / in equations of motion (37) to (40) by their equivalent expression given 
by (35) and (36) and averaging the resulting equations over the fast-changing angular 
variable /, the averaged equations will represent a system whose dynamics, in the first 
order of approximation, will be the same as the dynamical behavior of the original system 
at resonance, for a time interval e~ 1 / 2 t where t is the unit of time introduced in section 
2. To do this, let us differentiate equations (35) and (36) with respect to time and expand 
1/L 3 in powers of e 1 ! 2 as 



1 1 

~ Zjf 



(42) 



Equations (37) to (40) will now have the form 



D = -e 1/2 F u - eDF 12 + Q(e 3 / 2 ) , 



G 
9 



-eF 22 + 0(£ 3 / 2 ) , 
eF A2 + 0(e s / 2 ) , 



(43) 

(44) 
(45) 
(46) 



where 

Fn = ^H ext - AK L , (47) 
F 22 = ^-H ext — ATZq , (48) 

og 

F 3 2 = §i H e*t + Aft, , (49) 

^42 = T^ H ext + AU g , (50) 

and F 12 = dFn/dL are evaluated at (L , G , Lq 3 t + <p , (7) . In labeling these terms, I 
have followed the convention of Chicone et al. (1997b). 
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The averaged equations of the system at resonance are obtained by averaging equations 
(43) to (46) over the mean anomaly / . However, application of the averaging principle 
requires an appropriate averaging transformation that renders the system of (43) to (46) 
into a new system such that in the lowest order, it becomes exactly the first-order averaged 
system (Chicone et al. 1997b). To obtain this averaging transformation, one has to define 



X(G,<p,g,t) = Fu {G,-^t + <p,g,t) - F, 



li , 



(51) 



o 



where 




1-Kva I LO 



Fu(G,<p,g) 



Fu(G, —ot + <p,g,t) dt 



(52) 



o 



is the averaged value of Fn . Introducing A(G, <p, g, t) by 



d_ 
dt 



A(G,(p,g,t) = \(G,ip,g,t) 



(53) 



such that A = 0, one can use the averaging transformations D = D — e 1 ! 2 A(G, (p , g , t), 
G = G , <p = (p , g = g to write the equations (43) to (46) as 




) + 0(e 3 / 2 ) 



(54) 



G = -eF 22 + 0(e 3/2 ) , 



(55) 




6D 2 



(56) 



(57) 



Since equations (10) to (13) represent a perturbation problem of order e , I will keep the 
terms proportional to e in equations above and neglect the 0(e 3 ^ 2 ) terms. The averaged 
dynamics of the system will then be given by 
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D = -e 1/2 Fn - sDF 12 , (58) 
G = -eF 22 , (59) 

9 = eF 42 . (61) 



where F i2 , -F 2 2 , -F32 > and F 42 are calculated in the same fashion as F\\. 

6 FIRST ORDER AVERAGED DYNAMICS 

In the present work, only the dynamics of the first-order averaged system at reso- 
nance will be investigated. In this approximation, the first-order averaged equations have 
perturbation parameter e 1 / 2 and are given by 

D = -e^Fn , £ = (^) , 5 = ^=0, (62) 

where the terms with perturbation parameter s and higher have been neglected. 

Equation(62) is, in fact, a general representation of the first-order averaged system 
at resonance. In this equation, it is F\\ that contains information about the physical 
properties of the system. For the three-body system presented in this paper, F\\ is given 
by equation (47) where 

n L = --£-(i-e 2 )- 1/2 {RxC - R y v} . (63) 

In this equation 

C(L,G,l,g) = sin(t) + g) + e sing , (64) 

and 

T>(L,G,l,g) = cos(v + g) + e cos g . (65) 
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Since for the actual system of Sun- Jupiter-Saturn, numerical integrations have indi- 
cated a resonance lock with a near (2:1) commensurability (see section 3), one has to set 
m = 2 , n = 1 and u t = 1 , which is the angular frequency of Jupiter in dimensionless 
form, when one calculates Fn using integration (52). This integration requires H ext and 
71 l to be written in terms of the mean anomaly / . Appendices A and B contain details of 
these calculations. There, it has been shown that in the lowest order in eccentricity, TZl 
is proportional to e 2 . That is, 



where a ^ 1.625 represents the resonant value of the semimajor axis of Saturn (see Figure 
2). It is evident from this equation that the averaged value of TZl over the mean anomaly 
I will become zero. The first non- vanishing term in the averaged value of TZl appears as 
a term proportional to e 3 . Numerical value of this term is so small that its contribution 
in dynamics of the first-order averaged system becomes actually negligible (see Figure 2 ; 
the numerical value of the orbital eccentricity of Saturn at resonance is in average about 
0.022 with an amplitude of oscillation equal to 0.01). 

The appearance of the higher orders of eccentricity in equation (66) requires H ext to 
be expanded to higher orders in e. However, it turns out that even in expanding H ext 
to the second order in eccentricity, non-vanishing terms after averaging, are at most of 
order 10~ 4 . A comparison between these values and the order of magnitude of the term 
proportional to e in that expansion (i.e. 10 -2 ), indicates that the contributions of 0(e 2 ) 
terms are entirely negligible. Therefore, in the rest of these calculations, I will neglect 
the contribution of dynamical friction and will consider Fn ~ dH ext /dl. The external 
Hamiltonian H ext will also be expanded only to the first order in eccentricity. Appendix 
A shows that in that order, the only contribution of H ext in Fu appears as the term 
proportional to cos(2/ + g — 9j). Denoting this term by Hext, 




(66) 



Tlext — — 



ea 



cos(2Z + g - 6j) , 



(67) 
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where 



E 

h=0 



a&wr(§) 



1 + 



2/1 + 3 
/l + l 



3 /2/i + 5 



/i + 2 



4ag 



(68) 



From equation (52), Fn will now be equal to 



e a 



sin(2y? + g) , 



(69) 



and the first order averaged system at resonance will simply be obtained by substituting 
for F\\ in equation (62). That is, 



D 



a? 



e 1 ^ sin(2v? + g) , 



(70) 



and 



<p = — e 



1/2 [3D 

Us 



(71) 



where G and (7 are constants, and the tildes have been dropped for the sake of simplicity. 
From these two equations it follows that 



- 6eae ^ 

D — D cos(2y? + g) = , 

°o 



(72) 



and 



3eae , . 
<P — sin(2v? + g) = . 



a 



(73) 







Equation (73) is the equation of a mathematical pendulum. This equation can be 
attributed to a Hamiltonian in the form 



\pI + U{<p) 



(74) 



with a potential given by 



Uitp) 



Sea 

— j- cos(2(/? + g) + Const. 



2a 



(75) 
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Figure 7 shows the graph of U(<p) against <p . The periodic characteristic of U(<p) 
indicates that for all values of (p except for those corresponding to the maxima of U(ip) , 
the averaged system will fall into one of the potential wells and oscillates around a stable 
point forever. The resonance lock is established in this manner. For those values of (p 
where U (<p) becomes maximum, the system will be in an unstable equilibrium. A slight 
deviation from this instability will result in a resonance capture. 
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Figure 7. Graph of U (</?) against (p . The periodic nature of U (<p) guarantees resonance capture for essentially 
all initial relative positions of the two planets. 



Equation (72), too, shows a harmonic behavior (Figure 8) . It is important to remem- 
ber that D represents the deviation of L from its resonant value L . From the averaging 
principle presented in section 4, it is expected that D and L will have the same general 
behavior for time intervals of duration e~ 1 / 2 t . This is illustrated in Figure 8 for the system 
presented in this paper . 
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Figure 8. Graph of D (left) and the action variable L (right)agains time while the system is at resonance. The 
original system shows a periodic behavior (right) which in a time interval £ 1 t is in very good agreement 
with the first-order averaged system (left). The timescale is (10 4 T / /27T) years. 



7 SUMMARY AND CONCLUSIONS 



In an attempt to study analytically, the resonance capture phenomenon reported by 
Melita and Woolfson (1996), a restricted planar circular model of the Sun-Jupiter- Saturn 
system subject to dynamical friction with a freely rotating homogeneous interplanetary 
medium was studied. Numerical integrations of this system indicated a resonance lock 
with the same commensurability as reported by MW. The method of partial averaging 
was employed to study the dynamical evolution of the system while captured in resonance. 
By averaging the equations of motion over a fast-changing angular variable, these equations 
were reduced to the dynamical equations of a Hamiltonian system whose dynamics would 
be partially equivalent to the dynamical behavior of the main system. The application 
of this averaging method to the Sun- Jupiter-Saturn system where e = 0.001 , resulted 
in the first-order averaged system (70) and (71). Although this system is Hamiltonian 
and it is not expected to fully illustrate the dynamical evolution of the main dissipative 
system, but in a time interval s~ 1 / 2 T I /2n , it can present a reasonable approximation to the 
dynamics of the original system. The harmonic characteristic of the potential function of 
this Hamiltonian system guarantees the occurrence of resonance lock for all initial relative 
positions of the two planets. Also, the time evolution of the action variable of the first- 
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order averaged system illustrates an oscillatory behavior for the semimajor axis as well as 
the angular momentum of the main system over time scales of order e -1 / 2 . 

As mentioned above, the first-order averaged system at resonance is Hamiltonian. In 
general at this order, the dissipative effects of the non-gravitational perturbations appear 
as external torques in the equation of the mathematical pendulum presented by the dif- 
ferential equation of the angular variable ip. In the system presented in this paper, the 
averaged value of this torque in the lowest order in eccentricity, is proportional to e 3 and 
is entirely negligible when compared to the contribution of the gravitational attraction of 
the inner planet. The damping effect of this perturbation on the dynamical behavior of 
the original system can be illustrated by extending the calculations to the second-order 
partially averaged system. This will also allow for investigating the effects of the pertur- 
bation parameter e and the density of the interplanetary medium p Q on the time variation 
of the orbital elements of the outer planet and their resonant values (Haghighipour, in 
preparation) . 

The accreting force of the interplanetary medium was also neglected in this study. 
Although the density of the interplanetary medium was considered much higher than its 
present value, the masses of the planets were kept constant and equal to their present 
values. Numerical integrations performed by Melita and Woolfson (1996) have indicated 
that accretion will only change the results quantitatively. Therefore, one can apply the 
exact analysis presented in this paper, to the dynamical behavior of the system at resonance 
when accretion is also included. 

Numerical integrations were also carried out for different values of the perturbation 
parameter e and the density of the interplanetary medium p . The results indicated that 
when e is in the range 10~ 2 to 10~ 4 , a Sun- Jupiter-Saturn-like system will be captured 
in a near (2:1) resonance when the density of the interplanetary medium is taken to be 
equal to 16 times the mass of the inner planet spread uniformly in a sphere of radius 50 
au. However, for a real Sun- Jupiter-Saturn system in a less dense medium, the resonance 
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capture occurs for a shorter time interval and also moves to (3:1) commensurability. 

The ideas presented in this paper can be applied to any dynamical system at reso- 
nance. Of particular interest would be extension of calculations to include the gravitational 
tidal effects as another source of dissipation. This can be applicable to analyze the dynam- 
ical behavior of triple systems (Arzoumanian 1996) as well as newly discovered planetary 
systems (see for instance, Holland et al. 1998). 
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APPENDIX A : SIMPLIFICATION OF H, 



ext 



From equation (41) and in dimensionless form, H ext is given by 



H ext = - [r 2 - 2rcos(# - + l]' 



(Al) 



It is evident from this equation that in order to express H ext in terms of the mean anomaly 
/, one needs to write r and cos(6> — 6> 7 ) in terms of /. This is possible by writing r = 
G 2 (l + ecos-0) -1 and 9 = g + v, and expanding cos-0 and sinv as (Kovalevsky 1967) 



cos v = — e + 2 



i 2 00 

1 — e z 



and 



sint) = (1 - e 2 ) 1 / 2 ^ sin(jO [,/._, (je) - J j+1 (je)] 
In these equations, J is Bessel function of order j and is given by 



JM = E 



(-1)" 



^ v\(j + v)\ \2 



j+2v 



(A2) 



(A3) 



(A4) 



As mentioned in section 6, H ext is expanded to the first order in eccentricity. Substituting 
for cosv and sin-C in the expressions for r and cos(6> — 9 1 ) and considering terms of order 
e , one can therefore write 

r ~ a ( 1 — e cos /) , (A5) 



and 



cos(# - Oj) ~ cos(/ + # - 6>J + e cos(2/ + # - 0,) - cos(# - 0J ] . (A6) 



The external Hamiltonian H ext will therefore be given by 



H ext = - [l + a 2 - 2a cos(/ + g-9 I )\ 



! 1 2a 2 cosl + acos(2l + g - 6j) - 3acos(g - 9j) 
2 1 + a 2 -2acos(l + g -9j) 



(A7) 
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At resonance, the semimajor axis of the outer planet is almost constant with an 
average value ao — 1.625 . One can use this to simplify equation (A7) by expanding H ext 
in terms of (Iq 1 . The first term of this equation can be expanded using the relation 

\n - r 2 \ r> ^ Q Vr>/ 

where O is the angle between f\ and f*2 and P N (cos 0) is Legender polynomial of order 
N . In using equation (A8) in (A7), = / + g — Q I . Therefore the Legender polynomial 
P N (cos©) will produce terms proportional to cos N(l + g — Oj). Due to the harmonic 
nature of these terms, their averaged values obtained from equation (52) will all vanish. 
Therefore, the term in equation (A7) that is independent of the eccentricity e will have 
no contribution in F n . The remaining term in equation (A7) can also be expanded using 
the relation 

oo 

(1 - 2rcosa + r 2 )" A = £ (cosa) r q , |r| < 1 , (A9) 

q=0 

where (cos a) are Gegenbauer polynomials given by 

°i = £ ^i^-^y c ° s ^ - 2k) °] ■ (Ai °» 

In expanding the second term in equation (A7), A = 3/2, r = a^ 1 and a = I + g — Q 1 . 
From equation (36) where Lq = 2 , the only non-vanishing terms in integration (52) are 
the ones proportional to cos(2l + g — 9 I ) . Expansion (A10) produces terms proportional to 
cos[((/ — 2h) (I + g — Oj)] . Therefore after multiplying the numerator of the relevant term 
in equation (A7) by expansion (A9), the contribution of cos / appears in the terms where 
q = 2h + 1 , that of cos(2/ + g — 6j) in the terms where q = 2h and for cos(g — 6j) when 
q = 2h + 2 . The contributing terms in H ext (i.e., the terms with non- vanishing averaged 
values) can thus be written as in equation (67). 
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APPENDIX B : SIMPLIFICATION OF K L 



Substituting for R x and R y from equations (17) and (18) in (63), TZl can be written 

as 



n L = - 



BW 3 



(1 -e 2 )" 1/2 ln(l + Br 2 W 4 )\ a(l - e 2 ) (6 - u„) + ef sinv] . (Bl) 



The first step in simplifying TZl is to calculate the expression inside the bracket of 
equation (Bl). From equation (A5), the angular frequency of the medium uj^ = r -3 / 2 , 
and the angular frequency of the outer planet 9 = Gr~ 2 , are approximately equal to 



u>. 



~ a 3/2 ( 1 + ^ e cos / j , 



(B2) 



and 



6 ~ a" 3/2 (l + 2ecos/) , 



(B3) 



where G has been replaced by its equivalent expression from equation (27). Also from 
equation (A3) and to the first order in eccentricity, the radial momentum f = esinv/G 
can be written as 

r ~ ea~ 1/2 sin/ . (B4) 

Substituting these values in the expression inside the bracket of equation (Bl) and neglect- 
ing the terms of the second order and higher in eccentricity, 



a (1 — e 2 ) (9 — oj m ) + ef sinv 



1 



ea l / 2 cos/ 



(B5) 



It is now necessary to simplify the logarithmic term in equation (Bl). Numerical com- 
putations indicate that, when the system is at resonance (i.e. 10~ 4 < e < 10~ 2 ) , Br 2 W 4 
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is very small in comparison to 1 (Figure Bl). Therefore one can use the approximation 
ln(l + 5) ~ 5 for 5 « 1 to write this term as ln(l + Br 2 W A ) = Br 2 W 4 . 
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Figure Bl. The graph of Br 2 W 4 against time while the system is captured in resonance, a, e, 9, V and £ are 
respectively equal to 1.838046 , 0.0556 ,45 ,15 and 0.001 . This figure shows that at resonance, Br 2 W 4 IS 
almost constant with an order of magnitude at most equal to 10 — 2 . The timescale is (10 4 T J /27T) years. 



Replacing W 2 from equation (14) and using equations (A5) and (B2) to (B4), W can be 
written as 

W ~ a" 1 / 2 e(l - ^ cos 2 /) V2 . (B6) 

Finally from equations (B5) and (B6), TZl can be written as in equation (66) where in the 
lowest order, eccentricity appears as e 2 . 




